import plotly.graph_objects as go
import plotly.express as px
import numpy as np
import random
from scipy.spatial.transform import Rotation as R
class Filament:
def __init__(self, start, end):
self.start = start
self.end = end
def __repr__(self):
return "Start: {}, End: {}".format(self.start, self.end)
def BField(self, array_x, array_y, array_z, current):
"""
array_x, array_y, array_z must be ndarrays of equal shape with at least degree 3
I is a current value in amperes
"""
u0 = 1.25663706e-6
loop_array = np.array([self.start, self.end])
loop_diff = np.append(np.diff(loop_array, axis=0), [loop_array[0] - loop_array[-1]], axis=0)[:-1]
loop_shift = np.append(loop_array[1:], [loop_array[0]], axis=0)[:-1]
loop_array = loop_array[:-1]
point = np.array([array_x, array_y, array_z])
points = np.expand_dims(point, axis=1)
#### This tile function is fixed to (1, 2, 1, 1, 1) because there are only 2 points in a filament,
#### not a closed loop or multi segment
points = np.tile(points, (1, 2, 1, 1, 1))
points = points.T
# Vectors from points to endpoints
AP = points - loop_array
BP = points - loop_shift
# Field Math
r1 = np.sqrt((AP ** 2).sum(-1))[..., np.newaxis].T.squeeze().T
r2 = np.sqrt((BP ** 2).sum(-1))[..., np.newaxis].T.squeeze().T
Dot1 = np.multiply(AP, loop_diff).sum(-1)
Dot2 = np.multiply(BP, loop_diff).sum(-1)
cross = np.cross(loop_diff, AP)
CrossSqrd = (np.sqrt((cross ** 2).sum(-1))[..., np.newaxis]).squeeze() ** 2
top = (Dot1 / r1 - Dot2 / r2) * u0 * current
bottom = (CrossSqrd * 4 * np.pi)
factor = (top / bottom)
factor = factor[..., np.newaxis]
field = cross * factor
field = np.sum(field, axis=3)
field = field.T
return field[0], field[1], field[2]
def get_scatter3d(self, name, mode='lines+markers'):
line = go.Scatter3d(x=[self.start[0], self.end[0]],
y=[self.start[1], self.end[1]],
z=[self.start[2], self.end[2]],
name=name,
mode=mode)
return line
def get_cones(self, name, x_array, y_array, z_array, current):
U, V, W = U, V, W = newWire.BField(x_array, y_array, z_array, current)
cones = go.Cone(x=x_array.flatten(), y=y_array.flatten(), z=z_array.flatten(),
u=U.flatten(), v=V.flatten(), w=W.flatten(),
sizemode='scaled',
sizeref=1,
anchor='tail',
name=name)
return cones
# initialization function that takes in different parameters instead of a start, stop
def new_filament(start, length, deg_azimuth, deg_phi):
"""
Create a new filament at P_start with length oriented in the
direction specified by azimuth, phi
"""
start = np.array(start)
azi = deg_azimuth * np.pi / 180
phi = deg_phi * np.pi / 180
# Translate out of spherical coordinates
x,y,z = np.sin(azi), np.cos(azi), np.cos(phi)
dir_hat = np.array([x, y, z]) / np.linalg.norm(np.array([x, y, z]))
direction = length * dir_hat
endpoint = start + direction
return Filament(start, endpoint)
# Testing method a
new_filament([0,0,0], length=1.732, deg_azimuth=45, deg_phi=45)
# Testing method b
a = np.array([0,0,0])
b = np.array([1,1,1])
Filament(a, b)
a = np.array([0,0,0])
b = np.array([1,2,5])
newWire = Filament(a, b)
# Domain to draw the field over
X, Y, Z = np.meshgrid(np.linspace(-5,5,num=10),
np.linspace(-5,5,num=10),
np.linspace(-5,5,num=10))
line = newWire.get_scatter3d('NewLine')
cones = newWire.get_cones('Field', X, Y, Z, current=5)
fig = go.Figure(data=[cones, line], layout=dict(title=dict(text="B Field Plot")))
fig.show()
# Evaluate B Field at a point
i,j,k = 1,1,1
a,b,c = newWire.BField([[[i]]], [[[j]]], [[[k]]], current=5)
print(a.flatten(), b.flatten(), c.flatten())
# Declarations
TX_CURRENT = 1000
NUM_LINES = 10
# Domain to draw the field over
x,y,z = np.meshgrid(np.linspace(-10,10,num=20),
np.linspace(-10,10,num=20),
np.linspace(-10,10,num=20))
filaments = []
for i in range(NUM_LINES):
filaments.append(new_filament(start=[random.randint(-5,5),
random.randint(-5,5),
random.randint(-5,5)],
length=random.randint(1,10),
deg_azimuth=random.randint(0,359),
deg_phi=random.randint(0,180)))
fig = go.Figure(layout=dict(title=dict(text="Summed B Field Plot")))
Bx, By, Bz = np.zeros(x.shape), np.zeros(y.shape), np.zeros(z.shape)
for f in filaments:
fx, fy, fz = f.BField(x, y, z, current=TX_CURRENT)
fig.add_trace(f.get_scatter3d(str(random.randint(1,1000))))
Bx += fx
By += fy
Bz += fz
cones = go.Cone(x=x.flatten(), y=y.flatten(), z=z.flatten(),
u=Bx.flatten(), v=By.flatten(), w=Bz.flatten(),
sizemode='scaled',
sizeref=1,
anchor='tail',
name='Compiled Field')
fig.add_trace(cones)
fig.show()
## Visualize with streamlines
figstream = go.Figure(layout=dict(title=dict(text="Summed B Field Plot - Streamline")))
for f in filaments:
figstream.add_trace(f.get_scatter3d(str(random.randint(1,1000))))
# Initialize starting positions of streamlines
streamstarts = dict(x=[random.randint(-10,10) for i in range(100)],
y=[random.randint(-10,10) for i in range(100)],
z=[random.randint(-10,10) for i in range(100)])
streamlines = go.Streamtube(x=x.flatten(), y=y.flatten(), z=z.flatten(),
u=Bx.flatten(), v=By.flatten(), w=Bz.flatten(),
starts=streamstarts)
figstream.add_trace(streamlines)
figstream.show()
Discretize $T$ into N intervals $\frac{T}{N-1}$ wide as $dt$
N = 100
T = 1.0
dt = T/(N-1)
"""
Create a line using a random walk algorithm.
dt = arbitary t Interval
N = number of lines to generate
"""
dX = np.sqrt(dt) * np.random.randn(1, N)
dY = np.sqrt(dt) * np.random.randn(1, N)
dZ = np.sqrt(dt) * np.random.randn(1, N)
X = np.cumsum(dX, axis=1)
Y = np.cumsum(dY, axis=1)
Z = np.cumsum(dZ, axis=1)
points = np.vstack((X, Y, Z)).T
f_list = []
for i, p in enumerate(points[1:]):
f_list.append(Filament(points[i], p))
x,y,z = np.meshgrid(np.linspace(np.min(points[:,0]),np.max(points[:,0]),num=20),
np.linspace(np.min(points[:,1]),np.max(points[:,1]),num=20),
np.linspace(np.min(points[:,2]),np.max(points[:,2]),num=20))
fig = go.Figure(layout=dict(title=dict(text="Brownian Motion Simulation - BField")))
Bx, By, Bz = np.zeros(x.shape), np.zeros(y.shape), np.zeros(z.shape)
for fil in f_list:
fx, fy, fz = fil.BField(x, y, z, current=TX_CURRENT)
fig.add_trace(fil.get_scatter3d(str(random.randint(1,1000)), mode='lines'))
Bx += fx
By += fy
Bz += fz
cones = go.Cone(x=x.flatten(), y=y.flatten(), z=z.flatten(),
u=Bx.flatten(), v=By.flatten(), w=Bz.flatten(),
sizemode='scaled',
sizeref=1,
anchor='tail',
name='Compiled Field')
fig.add_trace(cones)
fig.show()
Define some survey geometry:
# Collar points of the boreholes
collar1 = np.array([550,550,0])
collar2 = np.array([700,850,0])
# e1 and e2 are the real space representations of the electrodes
e1 = np.array([600,600,-100])
e2 = np.array([700,800,-150])
# Add an extra borehole that we will "read" stations down
reading_collar = np.array([650,650,100])
reading_bottom = np.array([650,740,-200])
Generate some stations so we can profile the primary down the hole...
NUM_STATIONS = 40
# Get some stations between these points
v = reading_bottom - reading_collar
print("Borehole length: {}".format(np.linalg.norm(v)))
# Generate some stations
v_hat = v / np.linalg.norm(v)
stations = np.array([reading_collar + v_hat * x for x in np.linspace(0, np.linalg.norm(v), num=NUM_STATIONS)])
station_interval = np.linalg.norm(stations[1]-stations[0])
print("Length intervals: {}".format(station_interval))
station_names = np.linspace(0, np.linalg.norm(v), num=NUM_STATIONS)
stations
Set a hemispherical domain for each electrode, facing the other one to create a domain to iterate over for brownian movement
Pseudocode to generate a point from $P_1$ to $P_2$ where $P_1\neq P_2$, $P_{1z}\neq P_{2z}$ and $P_1$ and $P_2$ are position vectors in R3
We need to check the tolerance when walking to $P_2$. We will see if the generated point $P_n$ to $P_2$ is within tolerance. We will use a tolerance value where $|P_n - P_2| < Tolerance$ because generated floating point math is basically impossible to intersect exactly
Pseudocode:
function get_random_line:
create_a_ellipsoid_cone_domain
generate_a_point_in_domain
if point_in_tolerance (segment_length)
just append last point
return
else:
append_generated_point
recurse via get_random_line(generated_point)# RANGE LIMITS
AZIMUTH_RANGE = 85
PHI_RANGE = 85
PHI_RANGE *= np.pi/180
AZIMUTH_RANGE *= np.pi/180
def vector_ranges(start_p, end_p):
# Vector from e1 to e2
v = end_p - start_p
azi_rad = np.arctan2(v[0], v[1])
phi_rad= np.arccos(v[2] / np.linalg.norm(v))
# Because arctan resolves to values +-pi, we have to
# adjust for the negatives to get proper azimuths
if azi_rad < 0:
azi_rad += 2 * np.pi
azi_range = ((azi_rad + AZIMUTH_RANGE) % (2*np.pi),
(azi_rad - AZIMUTH_RANGE) % (2*np.pi))
# We will bind phi to the domain [0, 180]
# so that it doesn't wrap around
phi_range = (min(phi_rad + PHI_RANGE, np.pi),
max(phi_rad - PHI_RANGE, 0))
return azi_range, phi_range
a_r, p_r = vector_ranges(e1, e2)
print('Azimuth Range: {}'.format(np.array(a_r) * (180 / np.pi)))
print('Phi Range: {}'.format(np.array(p_r) * (180 / np.pi)))
Now we have the elliptic cone domain to call every step in the random walk, recall this function every step to regenerate the azimuth and phi ranges
def generate_random_points(start_p, target_p, segment_len):
azi_range, phi_range = vector_ranges(start_p, target_p)
azi_range = np.array(azi_range) * (180 / np.pi)
# Catch the wraparound for azimuth across North
if abs(azi_range[0] - azi_range[1]) > 180:
azi_range = [max(azi_range) - 360, min(azi_range)]
phi_range = np.array(phi_range) * (180 / np.pi)
newLine = new_filament(start=list(start_p),
length=segment_len,
deg_azimuth=random.uniform(azi_range[0], azi_range[1]),
deg_phi=random.uniform(phi_range[0], phi_range[1]))
distance = np.linalg.norm(newLine.end - target_p)
# Our tolerance distance will be the length of one segment
if distance < segment_len:
yield Filament(start_p, target_p)
else:
yield newLine
yield from generate_random_points(newLine.end, target_p, segment_len)
FILAMENT_LENGTHS = 20
# Test by generating paths between the two electrodes,
# with Filaments of length FILAMENT_LENGTHS
alistofwires = list(generate_random_points(e1, e2, FILAMENT_LENGTHS))
"Number of points in generated line: {}".format(len(alistofwires))
Plotting the created line
fig = go.Figure(layout=dict(title=dict(text="Random motion between two points")))
fig.add_trace(go.Scatter3d(x=[e1[0], e2[0]],
y=[e1[1], e2[1]],
z=[e1[2], e2[2]],
name="Line Between Electrodes",
mode='lines+markers'))
for f in alistofwires:
fig.add_trace(f.get_scatter3d(str(random.randint(1,1000)), mode='lines'))
fig.show()
TX_CURRENT = 2
# Define a new function that sums the fields in a list of filaments
def return_summed_BField(listofwires, x, y, z, current):
Bx, By, Bz = np.zeros(x.shape), np.zeros(y.shape), np.zeros(z.shape)
for fil in listofwires:
fx, fy, fz = fil.BField(x, y, z, current=current)
Bx += fx
By += fy
Bz += fz
return np.array([Bx, By, Bz])
# Get station coordinates
stnx = np.array(stations[:,0])
stny = np.array(stations[:,1])
stnz = np.array(stations[:,2])
BField_stations = [return_summed_BField(alistofwires,
np.array([[[sx]]]),
np.array([[[sy]]]),
np.array([[[sz]]]),
TX_CURRENT) \
for sx, sy, sz in zip(stnx, stny, stnz)]
# Cleanup
for i in range(len(BField_stations)):
BField_stations[i] = BField_stations[i].flatten()
BField_stations = np.array(BField_stations)
# Borehole vector
v
# For our example, dip is constant - straight line borehole
BH_dip = np.arctan2(v[2], np.linalg.norm(v[:2])) * -180/np.pi
BH_dip
# Borehole azimuth
BH_azimuth = np.arctan2(v[0], v[1]) * 180/np.pi
BH_azimuth
print('XYZ Primary BField at Stations')
BField_stations
Now we rotate from XYZ to AUV
# Rotate the theoretical values into the same frame of reference used with boreholes
rotatedBField = R.from_euler('Z', -90, degrees=True).apply(BField_stations)
# Rotate the theoretical values into the hole coordinate system
auvField = R.from_euler('YZ', [90 - BH_dip, BH_azimuth], degrees=True).apply(rotatedBField)
print('Rotated AUV Primary BField at Stations')
auvField
# Geometry Plot
x, y, z = [e1[0], e2[0]], [e1[1], e2[1]], [e1[2], e2[2]]
BH1 = [[e1[0], collar1[0]], [e1[1], collar1[1]], [e1[2], collar1[2]]]
BH2 = [[e2[0], collar2[0]], [e2[1], collar2[1]], [e2[2], collar2[2]]]
BHReading = [[reading_collar[0], reading_bottom[0]],
[reading_collar[1], reading_bottom[1]],
[reading_collar[2], reading_bottom[2]]]
fig = go.Figure()
fig.update_layout(title='Borehole MMR Profile Testing - Geometry')
fig.add_traces(go.Scatter3d(x=BH1[0], y=BH1[1], z=BH1[2], name='BH1'))
fig.add_traces(go.Scatter3d(x=BH2[0], y=BH2[1], z=BH2[2], name='BH2'))
fig.add_traces(go.Scatter3d(x=stations[:,0],
y=stations[:,1],
z=stations[:,2],
name='Survey BH',
hovertext=station_names,
mode='lines'))
for f in alistofwires:
fig.add_trace(f.get_scatter3d("Simulated Current Path", mode='lines'))
fig.show()
# AUV Profiling
Afig = go.Figure(layout=dict(title=dict(text="Primary Field Profile \"Z\"")))
Afig.add_trace(go.Scatter(x=station_names, y=auvField[:,0],
mode='lines+markers',
name='A Field'))
Afig.add_trace(go.Scatter(x=station_names, y=BField_stations[:,0],
mode='lines+markers',
name='Z - Unrotated BField'))
Afig.show()
Ufig = go.Figure(layout=dict(title=dict(text="Primary Field Profile \"X\"")))
Ufig.add_trace(go.Scatter(x=station_names, y=auvField[:,1],
mode='lines+markers',
name='U Field'))
Ufig.add_trace(go.Scatter(x=station_names, y=BField_stations[:,1],
mode='lines+markers',
name='X - Unrotated BField'))
Ufig.show()
Vfig = go.Figure(layout=dict(title=dict(text="Primary Field Profile \"Y\"")))
Vfig.add_trace(go.Scatter(x=station_names, y=auvField[:,2],
mode='lines+markers',
name='V Field'))
Vfig.add_trace(go.Scatter(x=station_names, y=BField_stations[:,2],
mode='lines+markers',
name='Y - Unrotated BField'))
Vfig.show()